import matplotlib.pyplot as plt
from skimage import io
import pandas as pd
import numpy as np
from numpy import random as rnd
Introduction:¶
Plaster models and synthetic modelling of plaster models is a way for geoscientists to get a better understanding of how processes in the subsurface work. When looking at seismic profiles, there might be a lack of understanding of how the structures in the seismic profiles were formed (Lecomte et al., 2015). A way to get a better understanding of how to interpret seismic profiles is through the use of synthetic seismic imaging being created from a plaster model. Here the geometrical evolution is known and there is a possibility to tweak the seismic modelling part to provide more realistic results as if the location was in the subsurface (Lindanger et al., 2004).
Figure 1: Simplification of data obtained.¶
For this term paper, fault lengths and dip-measurements of 5 plaster experiments have been acquired. The 5 plaster experiments were done with the same configuration to look at geometrical differences when an experiment is reproduced. A plaster model is a model created from liquid plaster. It is poured into a box that can either extend or compress the plaster. Ink is applied to the sides of the box, which are glass panes so that the compression/extension can be observed. When there is compression or extension, the ink will follow the motion of the plaster. The compression/extension takes place as the plaster starts to dry up, changing the material properties of the plaster, making brittle deformation possible. When the model has dried up faults can be seen in the model on both sides. These faults can then be measured to look at extension rate or compression rate, as well as fault characteristics under different regimes. These data can provide information regarding if the faults have similarities on both sides or if they had a different evolution. The plaster model can then be turned into synthetic seismic by using computer software. Figure 1 shows an illustrative explanation of how the data used is connected.
For this term paper regression and seismic image-analysis will be performed. The reasoning behind using the two data analysis methods is:
For regression, it would be interesting to explore the relationship between dip-measurements and fault lengths on both sides of the plaster models created. When doing this on both sides of the plaster models, it is possible to look at the correlation between the sides from a regression-analysis point of view.
The image analysis part of the term paper is an explorative way of comparing seismic sections. The reason why I want to do this approach is to make it easier for me to compare the different profiles, when there is different experiments being created from the same configuration. There are several steps to creating a synthetic seismic profile. The choices made during the processing will affect the final result, hence why it is important to being able to go from a result and comparing the result to the final imagery and approaches.
Beneath, in Figure 2, the seismic images processed from the same plaster model picture can be seen.
import matplotlib.image as mpimg
img_1 = mpimg.imread('./01_20Hz_45degDip.png')
img_2 = mpimg.imread('./02_40Hz_45degDip.png')
img_3 = mpimg.imread('./03_40Hz_30degDip.png')
img_4 = mpimg.imread('./04_40Hz_perfectIllumin.png')
fig, ax = plt.subplots(2,2, figsize =(40,20))
ax[0,0].imshow(img_1);
ax[1,0].imshow(img_2);
ax[0,1].imshow(img_3);
ax[1,1].imshow(img_4);
Figure 2: Synthetic seismic profiles created from a single plaster experiment.¶
At first, the data acquired needs to be imported into the .ipynb document. Beneath some explorative code can be seen. This code is included because it populates the data frame. This addition to the data frame gives a good overview of the length and dip differences in the dataset.
df = pd.read_excel("GEOV302_PlasterModel_Data.xlsx")
#Two columns needs to be assinged as numeric values in the dataframe. This has to do with how the data was set up in excel.
df['Maximum displacement, Side A (cm)'] = pd.to_numeric(df['Maximum displacement, Side A (cm)'],errors = 'coerce')
df['Maximum displacement, Side B (cm)'] = pd.to_numeric(df['Maximum displacement, Side B (cm)'],errors = 'coerce')
'''Calculating the dip-difference between fault initiation dip and dip of the fault at the end of the experiment. '''
dip_difference_a = []
dip_difference_b = []
for columns in df.itertuples(index = False):
if columns[3] > columns[2]:
diff_a = columns[3] - columns[2]
dip_difference_a.append(diff_a)
if columns[3] == columns[2]:
diff_a = 0
dip_difference_a.append(diff_a)
if columns[7] > columns[6]:
diff_b = columns[7] - columns[6]
dip_difference_b.append(diff_b)
if columns[7] == columns[6]:
diff_b = 0
dip_difference_b.append(diff_b)
if columns[7]<columns[6]:
diff_b = columns[6]-columns[7]
dip_difference_b.append(diff_b)
'''Calculating the difference in length of the same fault on both sides of the profile. '''
fault_length_difference = []
length_diff = 0;
for x in df.itertuples(index=True):
if x[5] < x[9]:
length_diff = x[9]-x[5]
fault_length_difference.append(length_diff)
if x[5] > x[9]:
length_diff = x[5]-x[9]
fault_length_difference.append(length_diff)
if x[5] == x[9]:
length_diff = 0;
fault_length_difference.append(length_diff)
'''Assigning the lists created earlier to the pandas dataframe.'''
df['Fault length difference (cm)'] = fault_length_difference
df['Difference in dip side A'] = dip_difference_a
df['Difference in dip side B'] = dip_difference_b
Table 1: A table outlining the data obtained from the plaster modelling in this experimental study on applying data-analysis applied with python on plaster model data.
df
Methods¶
Regression¶
To be able to perform a regression analysis using python, code has to be constructed. The code constructed can be seen in the cell below. This contains functions that have been created to make it possible to run a regression analysis on data columns of interest. At first, a function named regr_coeff was created. This calculates a regression coefficient by using the x- and y-values inputted. In the instance of applying this to the analysis of fault lengths/dip-measurements, the respective columns of the data of interest were put as x and y. Complimenting the function regr_coeff a function named yline_from_coeff was created. This function takes the intercept value and slope value created from regr_coeff combined with two new variables defined as xarray and xic to return the calculated regression line from the data. When these two functions where constructed, the code will address the outer limits of the interquartile range. This is done by calculating regression lines and ylines from the function several times over by using a random number generator. When these lines are produced, the confidence is set and the code produces an outer limit, set as 97,5% and 2,5%. Then the finalized regression plots for the maximum displacement of faults and dip angle at the end of the experiment are plotted against each other in a scatter plot, with the accompanying regression line and limits.
'''Defining the columns in the dataset, that I want to explore using regression analysis.'''
x = np.array(df['Maximum displacement, Side A (cm)'])
y = np.array(df['Dip angle at the end of experiment A'])
x2 = np.array(df['Maximum displacement, Side B (cm)'])
y2 = np.array(df['Dip angle at the end of experiment B'])
n=df.shape[0]
xic = x.min()
'''Defining a function that takes 3 inputs(x,y,xic). This function will do calculations to return the slope and interception
of the datapoints that is fed to the function.'''
def regr_coeff(x,y,xic):
xx = x*x
yy = y*y
xy = x*y
Sxx = n*(xx.mean() -x.mean() * x.mean())
Syy = n*(yy.mean() -y.mean() * y.mean())
Sxy = n*(xy.mean() -x.mean() * y.mean())
#From the values calculated above it is possible to determine the slope and intercept. These are then stored in the slope and icept variables.
slope = Sxy/Sxx
icept = np.mean(y) -slope*(np.mean(x) -xic)
return [icept,slope]
#Assigning slope and icept values as variables.
icept,slope = regr_coeff(x,y,xic)
xic2 = x.min()
icept2,slope2 = regr_coeff(x2,y2,xic2)
x_sorted = np.sort(x)
x_sorted2 = np.sort(x2)
'''Defining a function named yline_from_coeff that calculates the regression line from the datapoints given as the 4 inputs(icept,
slope,xarray,xic). icept and slope was calulated in the function defined above as the regr_coeff function.'''
def yline_from_coeff(icept,slope,xarray,xic):
yline = icept + slope*(xarray-xic)
return yline
#Calling the function:
yline = icept + slope*(x_sorted -xic)
yline2 = icept2 + slope2*(x_sorted2 -xic2)
#Defining a seed in random to get the same values each time this is computed.
rnd.seed(123)
x_sorted = np.sort(x)
x_sorted2 = np.sort(x2)
index = np.linspace(0,29,30).astype(int)
#Variable reps defines how many steps the bootstrap is to take.
reps = 5
#Defining a for-loop that runs from 0 to amount of repetitions(reps).
for i in range(0,reps):
#Defining a variable boot that is equal to random numbers generated in the dataset.
boot = rnd.choice(index, n, replace = True)
#Setting xr variable to x amount of boot.
xr = x[boot]
#Setting yr variable to y amount of boot.
yr = y[boot]
xr2 = x2[boot]
yr2 = y2[boot]
#Calculating slope and icept from the function created earlier.
(icept,slope) = regr_coeff(xr,yr,xic)
(icept2,slope2) = regr_coeff(xr2,yr2,xic2)
#Calculating the regressionline from the function 'yline_from_coeff' that was created earlier.
new_yline = yline_from_coeff(icept, slope, x_sorted,xic)
new_yline2 = yline_from_coeff(icept2, slope2, x_sorted2,xic2)
reps = 100
#defining a variable named arrayx, that is set to an array with the same size, but each value is set to zeros.
array_for_values = np.zeros([reps,n])
array_for_values2 = np.zeros([reps,n])
#Defining a for-loop that is set to run from 0 to the amount of reps.
for i in range(0, reps):
#Defining a variable boot that is equal to random numbers generated in the dataset.
boot = rnd.choice(index, n, replace = True)
#Setting xr variable to x amount of boot.
xr = x[boot]
#Setting yr variable to y amount of boot.
yr = y[boot]
#Calculating slope and icept from the function created earlier.
(icept,slope) = regr_coeff(xr,yr,xic)
#Calculating the regressionline from the function 'yline_from_coeff' that was created earlier.
new_yline = yline_from_coeff(icept, slope, x_sorted,xic)
#adding the values to arrayx from the function new_yline to fill it up with values calculated.
array_for_values[i,:] = new_yline[:]
xr2 = x2[boot]
#Setting yr2 variable to y2 amount of boot.
yr2 = y2[boot]
#Calculating slope2 and icept2 from the function created earlier.
(icept2,slope2) = regr_coeff(xr2,yr2,xic2)
#Calculating the regressionline from the function 'yline_from_coeff' that was created earlier.
new_yline2 = yline_from_coeff(icept2, slope2, x_sorted2,xic2)
#adding the values to arrayx2 from the function new_yline to fill it up with values calculated.
array_for_values2[i,:] = new_yline2[:]
#Setting quantile values(confidence interval end lines). To be plotted.
upper_boot = np.quantile(array_for_values, 0.975, axis =0)
lower_boot = np.quantile(array_for_values, 0.025, axis =0)
upper_boot2 = np.quantile(array_for_values2, 0.975, axis =0)
lower_boot2 = np.quantile(array_for_values2, 0.025, axis =0)
def regression_plot():
#Creating the initial figure.
fig = plt.figure(figsize =(12,6))
#Creating the subplots.
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#Plotting the scatter.
ax1.scatter(x,y)
ax2.scatter(x2,y2)
#Plotting the confidence interval end lines.
ax1.plot(x_sorted, upper_boot, color='pink',linestyle ='dashed',label='Upper confidence interval, 97.5%')
ax1.plot(x_sorted, lower_boot, color='pink',linestyle ='dashed', label = 'Lower confidence interval, 2.5%')
ax2.plot(x_sorted2, upper_boot2, color='pink',linestyle ='dashed',label='Upper confidence interval, 97.5%')
ax2.plot(x_sorted2, lower_boot2, color='pink',linestyle ='dashed', label = 'Lower confidence interval, 2.5%')
#Plotting the regression line.
ax1.plot(x_sorted,yline, marker='o', color = 'black', label='Regression line')
ax2.plot(x_sorted2,yline2, marker='o', color = 'black',label='Regression line')
#Adding title and labels to both subplots.
ax1.set_title('Fault length, plotted against final dip in side A')
ax1.set_xlabel('Fault length at end of experiment')
ax1.set_ylabel('Final dip at end of experiment')
ax2.set_title('Fault length, plotted against final dip in side B')
ax2.set_xlabel('Fault length at end of experiment')
ax2.set_ylabel('Final dip at end of experiment')
Image analysis¶
Two methods have been used to draw information from the pictures. One of the methods is then used to extract pixel values to quantify the data. The first method uses matplotlibs imread function to import the images. Then the red and blue bands are drawn from the picture. The blue band extracts troughs and the red band extract peaks. Then the red/blue-banded images are displayed next to each other. The reason for this is to get a better view of the images and their respective peaks/throughs. The second method uses skimage a scikit image processing package in python. Here the red and blue values representing peaks/troughs are extracted and plotted as histograms assorted into 9 bins. The intensity values of the respective images are extracted and saved as separate variables, which is then used to find a peak/trough ratio. The code for this part can be seen beneath.
seismic_image1 = plt.imread('./01_20Hz_45degDip.png')
seismic_image2 = plt.imread('./02_40Hz_45degDip.png')
seismic_image3 = plt.imread('./03_40Hz_30degDip.png')
seismic_image4 = plt.imread('./04_40Hz_perfectIllumin.png')
red1 = seismic_image1[:,:,0]
red2 = seismic_image2[:,:,0]
red3 = seismic_image3[:,:,0]
red4 = seismic_image4[:,:,0]
blue1 = seismic_image1[:,:,2]
blue2 = seismic_image2[:,:,2]
blue3 = seismic_image3[:,:,2]
blue4 = seismic_image4[:,:,2]
def color_band_images():
fig, ax = plt.subplots(4,2,sharex=True, sharey=True, figsize =(30,30))
ax[0,0].imshow(red1,cmap="bone")
ax[0,1].imshow(blue1, cmap="bone")
ax[0,0].title.set_text('20Hz 45degree dip, Red channel marked as white')
ax[0,1].title.set_text('20Hz 45degree dip, Blue channel marked as white')
ax[1,0].imshow(red2,cmap="bone")
ax[1,1].imshow(blue2, cmap="bone")
ax[1,0].title.set_text('40Hz 45degree dip, Red channel marked as white')
ax[1,1].title.set_text('40Hz 45degree dip, Blue channel marked as white')
ax[2,0].imshow(red3,cmap="bone")
ax[2,1].imshow(blue3, cmap="bone")
ax[2,0].title.set_text('40Hz 30degree dip, Red channel marked as white')
ax[2,1].title.set_text('40Hz 30degree dip, Blue channel marked as white')
ax[3,0].imshow(red4,cmap="bone")
ax[3,1].imshow(blue4, cmap="bone")
ax[3,0].title.set_text('40Hz Perfect Illumination, Red channel marked as white')
ax[3,1].title.set_text('40Hz Perfect Illumination, Blue channel marked as white')
fig.text(0.5, 0.04, 'Seismic image width', ha='center', size = 30)
fig.text(0.04, 0.5, 'Seismic image height', va='center', rotation='vertical',size= 30)
def hist_channels():
image1 = io.imread('./01_20Hz_45degDip.png')
image2 = io.imread('./02_40Hz_45degDip.png')
image3 = io.imread('./03_40Hz_30degDip.png')
image4 = io.imread('./04_40Hz_perfectIllumin.png')
fig, ax = plt.subplots(2,2, sharex=True, sharey=True, figsize =(20,10))
ax[0,0].hist(image1[:, :, 0].ravel(), bins = 9,color = 'red', alpha = 0.5)
ax[0,0].hist(image1[:, :, 2].ravel(), bins = 9, color = 'Blue', alpha = 0.5)
ax[0,0].legend(['Red_Channel', 'Blue_Channel'])
ax[0,0].title.set_text('20Hz 45degree dip')
ax[1,0].hist(image2[:, :, 0].ravel(), bins = 9, color = 'red', alpha = 0.5)
ax[1,0].hist(image2[:, :, 2].ravel(), bins = 9, color = 'Blue', alpha = 0.5)
ax[1,0].legend(['Red_Channel', 'Blue_Channel'])
ax[1,0].title.set_text('40Hz 45degree dip')
ax[0,1].hist(image3[:, :, 0].ravel(), bins = 9, color = 'red', alpha = 0.5)
ax[0,1].hist(image3[:, :, 2].ravel(), bins = 9, color = 'Blue', alpha = 0.5)
ax[0,1].legend(['Red_Channel', 'Blue_Channel'])
ax[0,1].title.set_text('40Hz 30degree dip')
ax[1,1].hist(image4[:, :, 0].ravel(), bins = 9, color = 'red', alpha = 0.5)
ax[1,1].hist(image4[:, :, 2].ravel(), bins = 9, color = 'Blue', alpha = 0.5)
ax[1,1].legend(['Red_Channel', 'Blue_Channel'])
ax[1,1].title.set_text('40Hz Perfect Illumination')
fig.text(0.5, 0.04, 'Intensity Value', ha='center', size = 30)
fig.text(0.04, 0.5, 'Count', va='center', rotation='vertical',size= 30)
red_count1 = np.sum(red1)
blue_count1 = np.sum(blue1)
ratio1 = np.sum(blue1)/np.sum(red1)
red_count2 = np.sum(red2)
blue_count2 = np.sum(blue2)
ratio2 = np.sum(blue2)/np.sum(red2)
red_count3 = np.sum(red3)
blue_count3 = np.sum(blue3)
ratio3 = np.sum(blue3)/np.sum(red3)
red_count4 = np.sum(red4)
blue_count4 = np.sum(blue4)
ratio4 = np.sum(blue4)/np.sum(red4)
name_list =['20Hz 45degree dip','40Hz 45degree dip','40Hz 30degree dip','40Hz Perfect Illumination']
red_count_list = [red_count1,red_count2,red_count3,red_count4]
blue_count_list = [blue_count1,blue_count2,blue_count3,blue_count4]
ratio_list = [ratio1,ratio2,ratio3,ratio4]
data = [red_count_list,blue_count_list,ratio_list]
ratio_df = pd.DataFrame()
ratio_df['Profile name'] = name_list
ratio_df['red count'] = red_count_list
ratio_df['blue count'] = blue_count_list
ratio_df['ratio'] = ratio_list
Results¶
From the code outlined in the methods part the plots and dataframe beneath has been produced.
Regression results:¶
The results from the regression analysis of the final dip at the end of the experiment and fault length at the end of the experiment can be seen in Figure 3. The first observation that is done is that the regression line for side A and side B has different linear trends, with side A dipping downwards and side B being steeper. Many of the scatter plots fall outside of the interquartile range that has been produced. There are some errors in the plots as well. Both regression plots have a scatter point that is situated at (0,0). This has to do with both sides having a fault that is not present on the other side. In regards to the two plot, it is important to observe that the respective plots does not show the same linear trend.
regression_plot()
Figure 3, regression plots of both sides with fault length plotted against final dip in respective side.¶
Image analysis results:¶
In Figure 4 the synthetically produced seismic images with their bands extracted have been plotted. Each row represents a single seismic image where the highlighted red band is shown as white on the left side. On the right side, the blue channel band representing a trough is highlighted as white. These images are inverse of each other. By looking at the different rows it shows that the frequenzy and degree of dip used to produce the seismic images creates different looking seismic imagery.
color_band_images();
Figure 4: Pictures where the color bands have been extracted and shown to highlight peak/trough visibility.¶
Figure 5 shows histograms. These histograms show intensity values of pixels plotted against pixel count. The histograms are made from the same model, which seismic images were produced from. They show the same general trend. Generally, for the first half of each graph, the dominant band is the colour blue. For the second half, the dominant band is the colour red. The things that differ in the graph is the ratio of blue-/red-bands that appears in the different bins. The purple colour shows where both the red and the blue channel is located. The count of the bins varies for the histograms, there are some bins that contain the majorities of the red and blue pixels in the seismic images, as well as the majority of the pixel count. In the second half, where red is dominating the majority of the count is either split between two bins or is located in one bin.
hist_channels();
Figure 5: Histograms of color intensity is plotted against pixel count of each image. Each histogram represents a different modelling approach.¶
From the red pixel count and blue pixel count values is extracted. A ratio between them is calculated. This ratio ranges from 0.9 to 0.93. The data can be seen in Table 2.
Table 2: Data frame of colourband-count.¶
ratio_df
Discussion¶
Regression: In the regression plots the trend of the regression line differs from the two sides. This is an indication that the relationship between fault length at the end of the experiment and fault dip at the end of the experiment differs from the two sides. This evidence can support that fault evolution on the two sides differ from each other, even though the faults created in the plaster is homogenous. This may be used as evidence to support that no fault is created perfectly planar.
Colour band images: This procedure was primarily to highlight the peaks and troughs in the seismic image. With this method, it is quite unclear if the picture displays areas where there is a very low peak or trough. This part was primarily not to highlight any form of data-analysis method, but more of a way to display data and familiarize with capabilities that python has in regards to showing images and image colour manipulation. If it was interesting to use these images further in the channel histograms this could be done by saving the pictures and using the images produced here to create the histograms. However, it is uncertain if the skimage python package would approve of this considering the bands had to be split even further, to represent white areas.
Channel histograms/Pixel ratios: The different histograms show a general trend. However, in some histograms, the red channel and blue channel count varies. It does not vary a lot, but it does show some difference. This part is of importance to highlight that even though there are only minor differences in pixel count of the images, there is an actual visible difference. A difference of 20Hz or a dip difference will pan out in a difference in the seismic imaging. This proves that the pixel ratios does catch differences in the seismic imagery. The count and ratios are quite similar. It is possible to see major differences in the histograms, but not in the ratios. This probably has to do with the very large number of pixels is counted. The count of pixels in some places display a difference of around 50 000, but this does not imply a great difference in the ratios, since the counts are such a large number.
Sources of error:¶
In the regression part of the term paper, a point of (0,0) appears in both plots. This has to do with a fault appearing on one side and not the other side of the plaster model. The reason for this is that the non-existent fault is registered in the data frame as 0 and not NaN. This does affect the regression analysis since the point has been used for calculating the regression. An important thing to consider is that the point appears in both plots and the correlation between them stays the same. It would be worse for the analysis if the (0,0) point only was present on one side, but not the other.
For the channel extraction histogram and following pixel ratios, there is a problem regarding if the analysis does not take into account the clarity of the reflectors. When doing just pixels it is simplified and pans down to just how many pixels are counted and not actually taking into account the geophysical parameters behind a seismic image. There might be even better ways to assess the reflectors, for example, clustering pixels to provide a more accurate result in terms of how geology is seen in a synthetic seismic image and how a computer sees it.
Final remarks¶
Unfortunately for me, data was not available through my supervisors. However, I managed to get a hold of data that shares some similarities with the work that I am going to perform in my master thesis. This for now is a framework that I will be able to continue when starting my actual master thesis works. I do think especially the image analysis part has some great potential when it comes to comparing different modelling approaches and choices during modelling. This term paper has provided me with some knowledge about image analysis and regression that I will bring with me further into the master thesis work, it also highlights the importance of how different approaches and how changes in modelling parameters will have a result for the final produced synthetic seismic.
References¶
I. Lecomte, P. L. Lavadera, I. Anell, S. J. Buckley, D. W. Schmid, M. Heeremans (2015). Ray-based seismic modelling of geologic models: understanding and analysing seismic images efficiently. Interpretation. Volume 3, No. 4, pp. SAC71-SAC89. doi: 10.1190/INT-2015-0061.1.
R. Lindanger, M. Øygaren, R. H. Gabrielsen, R. Mjelde, T. Randen, B. A. Tjøstheim. (2004). Analogue (plaster) modelling and synthetic seismic representation of hangingwall fault blocks above ramp-flat ramp extensional faults.R. Lindanger, M. Øygaren, R. H. Gabrielsen, R. Mjelde, T. Randen, B. A. Tjøstheim. (2004). Analogue (plaster) modelling and synthetic seismic representation of hangingwall fault blocks above ramp-flat ramp extensional faults.